Generic Chemistry
Generic chemical equilibrium and kinetic reactions following the standard approach used by reaction-transport codes such as PHREEQ and CrunchFlow, see eg (Steefel et al., 2015). This exploits timescale separation between "fast" (assumed instantaneous) chemical equilibrium reactions, and "slow" kinetic reactions or transport.
- The chemical system is represented by a small number of
totals
(or components) and an equal number of primary species concentrations, with secondary species concentrations calculated from the primary species via a set of equilibrium reactions. Primary species concentrations are determined by solving the set of algebraic equations given by the constraints on total concentrations.
- Kinetic reactions (with any species,
primaryorsecondaryas reactants) are then written withtotalsas products. - Bulk transport (eg ocean advection or eddy diffusivity) transports
totals. Molecular diffusivity (eg in a sediment) transportsprimaryorsecondaryspecies and accumulates fluxes intototals.
Reservoirs
PALEOaqchem.ReservoirsAq.ReactionConstraintReservoir — TypeReactionConstraintReservoirA primary species and (algebraic) constraint on a corresponding total or component.
The primary species concentration or amount is defined as a PALEO State Variable, which depending on the primary_variable parameter, may be:
Primary_conc: (mol m-3)Primary: (mol)Primary_pconc: -log 10 (concentration (mol kg-1))
The corresponding R_constraint_conc or R_constraint (mol) defining the algebraic constraint on the corresponding total (for use by the numerical solver) is defined as a PALEO Constraint Variable.
This ReactionConstraintReservoir would usually be used in combination with a ReactionReservoir that provides the required total component concentration or amount as an ODE variable (where as usual reaction source and sink fluxes are applied to the corresponding _sms variable). Depending on the constraint_variable parameter, the total component may be supplied as either a per-cell concentration or amount:
R_conc: (mol m-3)R: (mol)
Equilibrium reactions defining secondary species should add their contributions to the total to R_calc (mol). A primary species contribution R_calc += primary_total_stoich * Primary_conc * primary_volume is added to R_calc (where for the usual case parameter primary_total_stoich should be set to 1.0). Primary species contributions to other totals can be included by setting the primary_other_components parameter.
The numerical solver then solves for the primary species (and hence the secondary species concentrations) that (depending on the constraint_variable parameter) satisfy one of:
0 = R_constraint_conc = R_conc - R_calc/volume
0 = R_constraint = R - R_calcVolume conversions
The total species concentration R_conc and primary species concentration Primary_conc use (potentially different) volume conversions provided in volume and primary_volume respectively.
This allows for cases eg equilibrium partitioning between solute and solid phases by surface complexation, where R_conc refers to a cell total volume, and Primary_conc to a solute concentration.
Parameters
primary_total_stoich[Float64]=1.0,default_value=1.0,description="stoichiometric factor R_calc_conc += primary_total_stoich * Primary_conc"primary_other_components[Vector{String}]=String[],default_value=String[],description="contribution of primary species to other element or component total concentrations"primary_variable[String]="concentration",default_value="concentration",allowed_values=["concentration", "amount", "p_concentration"],description="units for primary variable"constraint_variable[String]="concentration",default_value="concentration",allowed_values=["concentration", "amount"],description="units for constraint variable"
Methods and Variables for default Parameters
do_constraintreservoir_primaryR_calc(mol),VT_ReactContributor,description="contributions to total R_calc_conc (NB: a total, not concentration, to generalize to multiphase eqb)"primary_volume–> volume (m3),VT_ReactDependency,description="cell volume (as used by Primary_conc)"Primary_conc(mol m-3),VT_ReactDependency,VF_State,description="concentration of primary species"
do_constraintreservoir_constraintR_calc(mol),VT_ReactTarget,description="contributions to total R_calc_conc (NB: a total, not concentration, to generalize to multiphase eqb)"R_constraint_conc(mol m-3),VT_ReactContributor,VF_Constraint,description="algebraic constraint on R_conc (= 0)"R_conc(mol m-3),VT_ReactDependency,description="total R_conc"volume(m3),VT_ReactDependency,description="cell volume (as used by total variable)"
PALEOaqchem.ReservoirsAq.ReactionImplicitReservoir — TypeReactionImplicitReservoirA primary species and corresponding total or component as an 'implicit' ODE variable.
This provides an implementation of the 'Direct Substitution Approach' to chemical speciation, where the total or component is a function of the primary species concentration.
The primary species concentration or amount is defined as a PALEO StateTotal Variable, which depending on the primary_variable parameter, may be:
Primary_conc: (mol m-3)Primary: (mol)Primary_pconc: -log 10 (concentration (mol kg-1))
The corresponding total component R_conc or R is defined as a PALEO Total Variable, which depending on the constraint_variable parameter, may be provided to the solver either as a per-cell concentration or amount:
R_conc = R_calc/volume: (mol m-3)R = R_calc: (mol)
Equilibrium reactions defining secondary species should add their contributions to the total to R_calc (mol). A primary species contribution R_calc += primary_total_stoich * Primary_conc * primary_volume is added to R_calc (where for the usual case parameter primary_total_stoich should be set to 1.0). Primary species contributions to other totals can be included by setting the primary_other_components parameter.
Source - sink fluxes eg kinetic reactions should be added to R_sms (mol yr-1) defined as a PALEO Deriv Variable.
Volume conversions
The total species concentration R_conc and primary species concentration Primary_conc use (potentially different) volume conversions provided in volume and primary_volume respectively.
This allows for cases eg equilibrium partitioning between solute and solid phases by surface complexation, where R_conc refers to a cell total volume, and Primary_conc to a solute concentration.
Parameters
primary_total_stoich[Float64]=1.0,default_value=1.0,description="stoichiometric factor R_calc_conc += primary_total_stoich * Primary_conc"primary_other_components[Vector{String}]=String[],default_value=String[],description="contribution of primary species to other element or component total concentrations"primary_variable[String]="concentration",default_value="concentration",allowed_values=["concentration", "amount", "pconcentration"],description="units for primary variable (specifies Primary\conc, Primary, Primary_pconc as StateTotal variable)"total_variable[String]="concentration",default_value="concentration",allowed_values=["concentration", "amount"],description="units for total variable (specifies R_conc, R as Total variable)"total[Bool]=false,default_value=false,description="true to calculate R_total = sum(R)"
Methods and Variables for default Parameters
do_implicitreservoir_primaryR_calc(mol),VT_ReactContributor,description="contributions to total R_calc_conc (NB: a total, not concentration, to generalize to multiphase eqb)"primary_volume–> volume (m3),VT_ReactDependency,description="cell volume (as used by Primary_conc)"Primary_conc(mol m-3),VT_ReactDependency,VF_StateTotal,description="concentration of primary species"
do_implicitreservoir_smsR_sms(mol yr-1),VT_ReactTarget,description="total or component R source - sink"R_conc_sms(mol m-3 yr-1),VT_ReactContributor,VF_Deriv,description="total or component R_conc source - sink"volume(m3),VT_ReactDependency,description="cell volume (as used by total variable)"
do_implicitreservoir_totalR_calc(mol),VT_ReactTarget,description="contributions to total R_calc_conc (NB: a total, not concentration, to generalize to multiphase eqb)"volume(m3),VT_ReactDependency,description="cell volume (as used by total variable)"R(mol),VT_ReactProperty,description="total or component R"R_conc(mol m-3),VT_ReactContributor,VF_Total,description="total or component R_conc"
PALEOaqchem.ReservoirsAq.ReactionAqConcSum — TypeReactionAqConcSumA sum of concentration variables (eg to get an element total)
Parameters
vars_to_add[Vector{String}]=["2\myvar", "myothervar", "-1\mythirdvar"],default_value=["2\myvar", "myothervar", "-1\mythirdvar"],description="vector of variable names to add, eg [2*myvar, myothervar, -1*mythirdvar]"add_to_sum_volume[Bool]=false,default_value=false,description="true to also add to a 'sum' Variable += 'sum_conc * volume"define_sum_volume[Bool]=false,default_value=false,description="only if 'add_to_sum_volume == true': true to also define the 'sum' Variable"
Methods and Variables for default Parameters
do_aqconcsumsum_conc(mol m-3),VT_ReactProperty,description="sum of specified variables"myvar(),VT_ReactDependency,description=""myothervar(),VT_ReactDependency,description=""mythirdvar(),VT_ReactDependency,description=""
Equilibrium reactions
PALEOaqchem.GenericReactions.ReactionAqEqb — TypeReactionAqEqbDefine a new equilibrium species N
N + a A <--> b B + c C
[N] = K_eqb'^K_power ([B]^b [C]^c) / ([A]^a)where to convert density units for K_eqb:
K_eqb' = K_eqb * rho_ref^K_density_powerThe first name in the Reactants list is the new species concentration N: other species concentrations in Reactants and Products lists must already be defined elsewhere in the model configuration.
The contribution of the new species to element totals or components is defined by the N_components vector, which may be empty eg to just calculate an Omega or a gas partial pressure etc.
Examples
Gas partial pressure from concentration
solubility_H2:
class: ReactionAqEqb
parameters:
Reactants: ["pH2"]
Products: ["H2_conc"]
K_eqb: 7.8e-1 # mol m-3 atm-1 at 298.15 K (Henry's law coefficent)
K_power: -1.0 # pH2 = H2_conc * K_eqb^-1
variable_attributes:
pH2%units: atmCompilation of Henry's law coefficients: https://www.henrys-law.org/ which is R. Sander: Compilation of Henry's law constants (version 5.0.0) for water as solvent, Atmos. Chem. Phys., 23, 10901-12440 (2023), doi:10.5194/acp-23-10901-2023
Unit conversions: 1 mol m^-3 Pa-1 = 1.01325e5 mol m-3 atm-1
Parameters
Reactants[Vector{String}]=["N\conc", "A\conc^2"],default_value=["N\conc", "A\conc^2"],description="concentrations or activities of new species followed by other reactants, write powers as X^2 etc"Products[Vector{String}]=["B\conc", "C\conc"],default_value=["B\conc", "C\conc"],description="concentrations or activities of products, write powers as X^2 etc"K_eqb[Float64]=0.00018629779999999998,default_value=0.00018629779999999998,description="equilibrium constant"K_density_power[Float64]=0.0,default_value=0.0,description="multiple K_eqb * rho_ref^K_density_power to convert units: 0.0 for K_eqb in mol m-3, 1.0 for K_eqb in mol kg-1, etc"K_power[Float64]=-1.0,default_value=-1.0,description="exponent of K_eqb"N_components[Vector{String}]=String[],default_value=String[],description="contribution of new species to element or component total eg '["2*TN_calc"]' to add2\*N\_conc\*volumetoTN\_calc(or empty vector to just define an Omega)"
Methods and Variables for default Parameters
do_aqeqbN_conc(mol m-3),VT_ReactProperty,description="aqueous concentration or activity"A_conc(mol m-3),VT_ReactDependency,description="aqueous concentration or activity"B_conc(mol m-3),VT_ReactDependency,description="aqueous concentration or activity"C_conc(mol m-3),VT_ReactDependency,description="aqueous concentration or activity"
Kinetic reactions
PALEOaqchem.GenericReactions.ReactionAqKinetic — TypeReactionAqKineticDefine a kinetic reaction with rate dependent on concentrations
a A + b B --> c C + d DRate (default) is:
R = K * [A] * [B]where this can be modified to different functional form by defining a vector of Rate_functions to apply to each concentration variable.
Parameters Reactants and Products should be the vectors of stoichiometry * <name> of (total) species to accumulate fluxes into <name>_sms variables.
Parameter Reactant_concs should be an empty vector to take default concentration variable names from Reactants, or a Vector of names to specify concentration species names explicitly (required when eg where Reactants refers to totals which are partioned into multiple species).
Parameters
Reactants[Vector{String}]=["A", "2\B"],default_value=["A", "2\B"],description="reactant species"Products[Vector{String}]=["2\C", "D"],default_value=["2\C", "D"],description="product species"Reactant_concs[Vector{String}]=String[],default_value=String[],description="names of concentration variables to calculate rate eg '[`"A_conc"]' etc, empty vector to used defaults from 'Reactants' eg 'A_conc', 'B_conc' ..."Rate_functions[Vector{String}]=String[],default_value=String[],allowed_values=["linear", "sqrt", "monod"],description="functional form for rate function of each concentration (empty vector for default 'linear')"K[Float64]=NaN,default_value=NaN,description="rate constant"K_lim[Float64]=NaN (mol m-3),default_value=NaN,description="limiting concentration for 'monod' rate function"
Methods and Variables for default Parameters
do_aqkineticredox_A_B_C_D(mol yr-1),VT_ReactProperty,description="rate variable"A_conc(mol m-3),VT_ReactDependency,description="aqueous concentration or activity"B_conc(mol m-3),VT_ReactDependency,description="aqueous concentration or activity"volume(m3),VT_ReactDependency,description="cell solute volume"
RateStoich_redox_A_B_C_Dredox_A_B_C_D(mol yr-1),VT_ReactDependency,description="rate variable"- [
A_sms] (mol yr-1),VT_ReactContributor,description="generated by RateStoich rate=redox_A_B_C_D" - [
B_sms] (mol yr-1),VT_ReactContributor,description="generated by RateStoich rate=redox_A_B_C_D" - [
C_sms] (mol yr-1),VT_ReactContributor,description="generated by RateStoich rate=redox_A_B_C_D" - [
D_sms] (mol yr-1),VT_ReactContributor,description="generated by RateStoich rate=redox_A_B_C_D"
totalsredox_A_B_C_D(mol yr-1),VT_ReactDependency,description="rate variable"redox_A_B_C_D_total(mol yr-1),VT_ReactProperty,description="total rate variable"
Precipitation-dissolution reactions
PALEOaqchem.GenericReactions.ReactionAqPrecipDissol — TypeReactionAqPrecipDissolDefine a precipitation and dissolution reaction for solid S
a A + b B <--> s S + d DRate for the precipitation and dissolution reactions are:
R_precip = K_precip * (Ω - 1) (Ω > 1)
R_dissol = K_dissol * S_conc * (1 - Ω) (Ω < 1)Parameters Reactants and Products should be the vectors of stoichiometry * name of (total) species to accumulate fluxes into _sms variables.
Solid_conc should be the name of the concentration variable for S, or an empty string to use default 'S_conc'.
Parameters
Reactants[Vector{String}]=["A", "2\B"],default_value=["A", "2\B"],description="reactant species"Products[Vector{String}]=["S", "0.5\D"],default_value=["S", "0.5\D"],description="product species (solid S first)"Solid_conc[String]="",default_value="",description="name of solid S concentration variable (empty string to default to 'S_conc')"K_precip[Float64]=0.0 (mol m-3 yr-1),default_value=0.0,description="rate constant for precipitation reaction"K_dissol[Float64]=0.0 (yr-1),default_value=0.0,description="rate constant for dissolution reaction"dissol_rolloff_conc[Float64]=0.0 (mol m-3),default_value=0.0,description="limiting concentration below which dissolution rate is rolled off to zero as Solid_conc^2"
Methods and Variables for default Parameters
do_aqprecipdissolprecipdissol_A_B_S_D(mol yr-1),VT_ReactProperty,description="rate variable"S_conc(mol m-3),VT_ReactDependency,description="solid concentration or activity"Omega(),VT_ReactDependency,description="saturation state"volume(m3),VT_ReactDependency,description="cell solid phase volume"
RateStoich_precipdissol_A_B_S_Dprecipdissol_A_B_S_D(mol yr-1),VT_ReactDependency,description="rate variable"- [
A_sms] (mol yr-1),VT_ReactContributor,description="generated by RateStoich rate=precipdissol_A_B_S_D" - [
B_sms] (mol yr-1),VT_ReactContributor,description="generated by RateStoich rate=precipdissol_A_B_S_D" - [
S_sms] (mol yr-1),VT_ReactContributor,description="generated by RateStoich rate=precipdissol_A_B_S_D" - [
D_sms] (mol yr-1),VT_ReactContributor,description="generated by RateStoich rate=precipdissol_A_B_S_D"
totalsprecipdissol_A_B_S_D(mol yr-1),VT_ReactDependency,description="rate variable"precipdissol_A_B_S_D_total(mol yr-1),VT_ReactProperty,description="total rate variable"
Particulate fluxes
PALEOaqchem.Particle.ReactionParticleDecay — TypeReactionParticleDecayDecay (eg organic matter to remineralization) at rate 1.0/decay_timescale of eg an organic matter dissolved/particulate phase in Reservoir Particle, to decayflux. Particle may have isotope_type. The Reservoir Particle may have the :vsink attribute set to represent a marine sinking particulate phase.
Parameters
decay_timescale[Float64]=0.5 (yr),default_value=0.5,description="particle decay timescale"decay_threshold[Float64]=-Inf (mol m-3),default_value=-Inf,description="particle decay concentration below which decay stops"show_decayrate[Bool]=false,default_value=false,description="true to create a 'decayrate' variable to record decay rate"field_data[DataType]=PALEOboxes.ScalarData,default_value=PALEOboxes.ScalarData,allowed_values=Type[PALEOboxes.ScalarData, PALEOboxes.IsotopeLinear],description="disable / enable isotopes and specify isotope type"
Methods and Variables
do_particle_decayParticle(mol),VT_ReactDependency,description="Particle amount"Particle_sms(mol yr-1),VT_ReactContributor,description="Particle source-sink"decayflux(mol yr-1),VT_ReactContributor,description="Particle decay flux"
PALEOaqchem.Particle.ReactionFluxToComponents — TypeReactionFluxToComponentsDistribute a single input_flux (eg an organic matter flux) to output_fluxnames components with stoichiometry output_fluxstoich. input_flux may have an isotope type (set by field_data) in which case the isotopic composition will be sent to (usually one) output_fluxname with ::Isotope suffix.
Parameters
outputflux_prefix[String]="",default_value="",description="Prefix for output flux names"outputflux_names[Vector{String}]=["Corg", "N", "P"],default_value=["Corg", "N", "P"],description="Suffixes for output flux names. Use ::Isotope suffix to identify a flux with 'isotope_type'"outputflux_stoich[Vector{Float64}]=[106.0, 16.0, 1.0],default_value=[106.0, 16.0, 1.0],description="Stoichiometry for output fluxes relative to input flux"allow_unlinked_outputflux[Bool]=false,default_value=false,description="true to allow (and ignore) unlinked outputflux Variables, false to error if any outputflux Variable is not linked"field_data[DataType]=PALEOboxes.ScalarData,default_value=PALEOboxes.ScalarData,allowed_values=Type[PALEOboxes.ScalarData, PALEOboxes.IsotopeLinear],description="disable / enable isotopes and specify isotope type"
Methods and Variables for default Parameters
do_flux_to_componentsinputflux(mol yr-1),VT_ReactTarget,description="input flux"Corg(mol yr-1),VT_ReactContributor,description="flux Corg"N(mol yr-1),VT_ReactContributor,description="flux N"P(mol yr-1),VT_ReactContributor,description="flux P"
Co-precipitation
PALEOaqchem.CoPrecip.ReactionPACoPrecip — TypeReactionPACoPrecipCo-precipitation of P (eg iron-bound phosphorus) with A (eg Fe oxide) formation
P -> P=Aat a rate gamma * A_formation_rate_<n>, with limitation at low P concentration
P_components defines solution totals or components that should be modified when P is consumed (eg ["-1*P"] to remove "P" from solution)
Parameters
A_rate_stoich_factors[Vector{Float64}]=[1.0],default_value=[1.0],description="stoichiometry factor to multiply each A formation rate variable to convert to mol A"gamma[Float64]=0.15 (mol/mol),default_value=0.15,description="P:A molar ratio"P_limit[Float64]=1.0e-6 (mol m-3),default_value=1.0e-6,description="limiting P concentration below which co-precipitation is inhibited"P_components[Vector{String}]=["-1\P", "TAlk"],default_value=["-1\P", "TAlk"],description="totals or components that should be modified when P is consumed from solution"
Methods and Variables for default Parameters
do_PA_copreciprate_PA_coprecip(mol P yr-1),VT_ReactProperty,description="rate of P co-precipitation"P_conc(mol m-3),VT_ReactDependency,description="P concentration"A_formation_rate_1(mol m-3 yr-1),VT_ReactDependency,description="substance A formation rate"
RateStoich_rate_PA_copreciprate_PA_coprecip(mol P yr-1),VT_ReactDependency,description="rate of P co-precipitation"- [
P_sms] (mol yr-1),VT_ReactContributor,description="generated by RateStoich rate=rate_PA_coprecip" - [
TAlk_sms] (mol yr-1),VT_ReactContributor,description="generated by RateStoich rate=rate_PA_coprecip" - [
PA_sms] (mol yr-1),VT_ReactContributor,description="generated by RateStoich rate=rate_PA_coprecip"
totalsrate_PA_coprecip(mol P yr-1),VT_ReactDependency,description="rate of P co-precipitation"rate_PA_coprecip_total(mol P yr-1),VT_ReactProperty,description="total rate of P co-precipitation"
PALEOaqchem.CoPrecip.ReactionPARelease — TypeReactionPAReleaseRelease of P (eg iron-bound phosphorus) with A (eg Fe oxide) destruction
P=A -> Pat a rate Prelease = theta * A_destruction_rate_<n>, where theta = PA_conc / A_conc
P_components defines totals or components that should be modified when P is released (eg ["P"] to return P to solution).
Parameters
A_rate_stoich_factors[Vector{Float64}]=[1.0],default_value=[1.0],description="stoichiometry factor to multiply each A destruction rate variable to convert to mol A"P_components[Vector{String}]=["P", "-1\TAlk"],default_value=["P", "-1\TAlk"],description="totals or components that should be modified when P is released"
Methods and Variables for default Parameters
do_PA_releaserate_PA_release(mol P yr-1),VT_ReactProperty,description="rate of coprecipitated P dissolution"PA_conc(mol m-3),VT_ReactDependency,description="adsorbed P concentration"A_conc(mol m-3),VT_ReactDependency,description="adsorbant concentration"PA_theta(mol/mol),VT_ReactProperty,description="P:A molar ratio of adsorbed of coprecipitated P"A_destruction_rate_1(mol m-3 yr-1),VT_ReactDependency,description="substance A destruction rate"
RateStoich_rate_PA_releaserate_PA_release(mol P yr-1),VT_ReactDependency,description="rate of coprecipitated P dissolution"- [
PA_sms] (mol yr-1),VT_ReactContributor,description="generated by RateStoich rate=rate_PA_release" - [
P_sms] (mol yr-1),VT_ReactContributor,description="generated by RateStoich rate=rate_PA_release" - [
TAlk_sms] (mol yr-1),VT_ReactContributor,description="generated by RateStoich rate=rate_PA_release"
totalsrate_PA_release(mol P yr-1),VT_ReactDependency,description="rate of coprecipitated P dissolution"rate_PA_release_total(mol P yr-1),VT_ReactProperty,description="total rate of coprecipitated P dissolution"